#include <cmath>
#include <iostream>
#include <cstdlib>
using namespace std;
#include "RateLaw_Kleijn1989.h"

/*****************************************************************************/
double InvReactionModel(double w, double P, double T)
{
	double w_max = ReactionModel(1.0,P,T);
	if (w<0.0||w>w_max)		{cerr << "ERROR@InvReactionModel: Invalid input. w(mole/m^2*s) = " << w << "\n"; system("pause");exit(0);}
	if (P<0.0)				{cerr << "ERROR@InvReactionModel: Invalid input. P(Pa) = " << P << "\n"; system("pause");exit(0);}
	if (T<0.0||T>1500)		{cerr << "ERROR@InvReactionModel: Invalid input. T(K) = " << T << "\n";system("pause");exit(0);}

	// w: layer/s
	double Xs = 0.001;
	double w_temp;
	
	w_temp = ReactionModel(Xs,P,T);
	while(w_temp<w)
	{
		Xs+=0.001;
		w_temp = ReactionModel(Xs,P,T);
	}
	if(Xs>1.0)
	{
		cerr << "ERROR@InvReactionModel: W is too large. w(layer/s) = " << w << "\n";
		system("pause");exit(0);
	}

	return Xs;
}
/*****************************************************************************/
double ReactionModel(double Xs, double P, double T)
// INPUTS:
//   Xs: molar fraction
//   P : pressure, Pa
//   T : Temperature, K
// OUTPUT
//   Reaction rate, layer/s
{
	if (Xs < 0.0 || Xs>1.0) {cerr << "ERROR@ReactionModel: Invalid input. Xs = " << Xs << "\n";system("pause");exit(0);}
	if (P<0.0)				{cerr << "ERROR@ReactionModel: Invalid input. P(Pa) = " << P << "\n";system("pause");exit(0);}
	if (T<0.0||T>1500)		{cerr << "ERROR@ReactionModel: Invalid input. T(K) = " << T << "\n";system("pause");exit(0);}

	return (5.7423e4)*(1.6*10000.0*exp(-18500/T)*(P*Xs))/(1+0.19*sqrt((1-Xs)*P)+0.70*(P*Xs)); // layer/s
}